############################################
############################################
# previous: 7-VAP_FutSum_Thresh_Plot########
############################################

###### Summary ############################################################################################
#This script calculates the cost distance from known occurrences for each species based on habitat suitability for their
#modeling unit's current, future median, Q10, and Q90 rasters. Cost distance is calculated assuming that ability to spread 
#from known occurrences into surrounding habitat is limited by the suitability of each grid cell for the species, 
#i.e. suitable cells are easier and more likely to function as corridors for species dispersal than unsuitable cells.
#The species is assumed to be more likely to be present in suitable habitat that is connected to known occurrences, 
#even if it hasn't yet been observed there, but unlikely to be present in disconnected habitat patches or moderately 
#suitable habitat far away from known occurrences. Cost distance is calculated based on un-thresholded suitability, 
#i.e. even cells with a suitability value below 'presence' level are considered as dispersal corridors of high resistance.
#Habitat suitability models are cut to a cost distance value up to x (='cost distance' kilometers rather than linear kilometers).
#models are then further cut to within X linear distance km of known occurrences. Resulting rasters are thresholded and
#saved as continuous (habitat suitability 0-1) and binary (presence absence 1/0) rasters. Future (median, Q10, Q90) rasters are 
#furthermore restricted to within a Xkm linear distance buffer of final current, cost-distance cut presence rasters to represent 
#the assumption that species cannot expand their ranges by more than Xkm compared to current distributions by 2050
#even if suitable habitat expands further because of population dispersal and establishment limitations. All rasters created in 
#this analysis process are saved and .png. maps of outputs created. Note that output files are created for each species, not for 
#each modelling unit. While species in some cases had to be combined into modeling units to allow for statistical rigor of models,
#the final outputs created here are based on the suitable habitat reachable by each individual species based on the location of 
#their known occurrences and the habitat suitability of their overall modeling unit. The cost-distance cut outputs represent 
#likely area of occurrence in contrast to the simply thresholded but un-cut previous outputs, which simply show where the species 
#'could' occur rather than where it is 'known to' (actual occurrences) or 'likely to' (reachable areas of sufficient suitability) occur.
#
#divide likely to occur habitat into core (e.g. x>=0.5) and non-core (e.g. threshold<x<0.5) habitat?
###########################################################################################################

############################################
# next: 9-VAP_Hotspots #####################
############################################

#################################################################################################################
##### create cost distance and restricted distribution models ###################################################
#################################################################################################################

#cut modeling units back into individual species to allow for visualization of each separate species 
#even if it does not have enough records to create its own model, then merge into species groups: closely 
#related species with similar venom that are unlikely to have additive abundances: i.e. they probably don't 
#overlap and if they do they probably form transition zones or patchworks from one species to the other 
#rather than having cumulative abundances. This avoids overestimating diversity and abundance - hotspots 
#can then be calculated by adding distinct species groups together (e.g. carpet vipers and spitting cobras).

#load packages:
library(R.utils)
library(raster)
library(maptools)
library(rgdal)
library(sp)
library(stringr)
library(gdistance)
library(Rcpp)
library(sf)

#Main working dir:
MyDir<-paste("/home/jc217070/Maxent",sep="")  #define working directory on HPC
OccDir<-file.path(MyDir,"Maxent_SWDs/spp_occs")
SWDDir<-file.path(MyDir,"Maxent_SWDs/spp_swds")
BGDir<-file.path(MyDir,"Maxent_SWDs/backgrounds")
FullSDMDir<-paste(MyDir,"/Outputs1_Full",sep="")
OutDirCrossval<-paste(MyDir,"/Outputs2_Crossval",sep="")
OutDirFin<-file.path(MyDir,"Outputs3_Final")
PlotDir<-paste(MyDir,"/Maxent_Data/Plots/Plots_WHOSnakes_2020_09_01",sep="")
IndSppOutDir<-paste(OutDirFin, "/Final_Individual_Spp", sep="")

#set project name:
projectName<-"WHOSnakes"

#define selection criteria of interest for further analyses:
crit<-c("permutation.importance")   #c("permutation.importance","contribution","Test.gain.without","Test.gain.with.only","AUC.without","AUC.with.only")
critThresh<-c(1)                     #c(1,1,1,10,1,75)

#define maximum cost distance allowed from known occurrences (300km=300000M)

x1 <- 250000
x2 <- 500000

#read species reference list:
DatSum<-read.table(file=file.path(paste(MyDir, "/Maxent_LUTables/",projectName,".csv",sep="")), header=TRUE, sep=",") #load data summary for all species
#DatSum<-read.table(file=file.path(paste(MyDir, "/Maxent_LUTables/",projectName,".csv",sep="")), header=TRUE, sep=",") #load data summary for all species
spp<-DatSum$gen_sp_subtax

mapspp<-unique(DatSum$ModUnit)
mapspp<-paste(gsub(" ","_",mapspp))

mapsppB<-unique(DatSum$ModUnit)
# mapsppB<-unique(DatSum$ModUnit[DatSum$MergeShort %in% c("Naja naja complex",
#                                                         "Daboia russelii group",
#                                                         "Bungarus caeruleus complex",
#                                                         "Echis carinatus",
#                                                         "Echis coloratus super-complex")])
#mapsppB<- (c("Bothrops asper","Bothrops atrox"))

sppB<-DatSum$gen_sp_subtax[DatSum$ModUnit %in% mapsppB]
mapsppB<-paste(gsub(" ","_",mapsppB))

#define GDA94 gcs projection that I want to project rasters to
proj4_WGS84 <- CRS("+init=epsg:4326")
#make 1km snake raster:

bg<-raster(paste(PlotDir,"/toprug_1km/w001000.adf",sep=""), crs="+init=epsg:4326")
GLOBraster<-bg
values(GLOBraster)<-ifelse(values(GLOBraster)!=-Inf,0,-Inf) #raster showing all land as 0, ocean as NA
GLOBrasterB<-GLOBraster
values(GLOBrasterB)<-ifelse(values(GLOBrasterB)==0,1,values(GLOBrasterB)) #raster showing all land as 1, ocean as NA
GLOBrasterC<-GLOBraster
GLOBrasterC[GLOBrasterC == 0] <- NA

worldmap <- st_read(paste(PlotDir,"/worldmap/WHO_countries.shp",sep=""))
WHOpolies <- st_read(paste(MyDir,"/WHOSnakesShapes/SSN_CTY_CAT_BRG_FEA_IN_Mar2023.shp",sep=""))

#outcats<-c("Current","2050_Median","2050_Q90","2050_Q10","2090_Median","2090_Q90","2090_Q10")
outcats<-c("2050_Median","2090_Median","2050_Q90","2050_Q10","2090_Q90","2090_Q10")

rasterOptions()

# #################################################################################
# #make temporary temp dir for rasters in this job and set from default to this dir
# dir.create(paste(MyDir,"/tmpCutCD",sep=""),showWarnings=FALSE)
# rasterOptions(tmpdir=paste(MyDir,"/tmpCutCD",sep=""))
# #################################################################################

for (cat in outcats){
  
  ###################################
  # general mapspp info #############
  ###################################
  
  for (mapsp in mapsppB){

    ########################################
    #delete any temp files older than 1h
    removeTmpFiles(h=0.1)
    ########################################
    
    is<-which(DatSum$ModUnit==gsub("_"," ",mapsp))
    myspp<-DatSum$gen_sp_subtax[is]
    mapspnum<-which(mapspp==mapsp)
    
    
    
    #only run for this mapsp if the last species included in the unit hasn't been processed previously:
    if(!file.exists(paste(IndSppOutDir,"/",myspp[length(myspp)],"_",cat,"_CropThreshCDCutBin.asc",sep="")) &
       !file.exists(paste(IndSppOutDir,"/",myspp[length(myspp)],"_",cat,"_CropThreshCDCutBin.asc.gz",sep=""))){
    
    
    
    
    SpOutDirFin<-paste(OutDirFin, "/", crit,"/ModUnit_",mapsp, sep="")
    SpOutFile<-paste(SpOutDirFin, "/ModUnit_", mapsp, "_", cat, "Crop.asc", sep="")
    if(cat != "Current"){
      SpOutFile<-paste(SpOutDirFin, "/future/ModUnit_", mapsp, "_", cat, "Crop.asc", sep="")
    }
    
    
    if (file.exists(paste(SpOutFile,".gz", sep="")) & !(file.exists(paste(SpOutFile)))){
      gunzip(filename=paste(SpOutFile,".gz",sep=""),
             destname=paste(SpOutFile),
             remove=FALSE, overwrite=TRUE)
    }
    
    if (file.exists(paste(SpOutFile))){
      
      # read in the maxent results thresholds
      results <- read.csv(paste(SpOutDirFin, "/maxentResults.csv", sep=""))
      TrainAUC<-results$Training.AUC[1]
      resultsCV <- read.csv(paste(OutDirCrossval, "/", crit, "/ModUnit_", mapsp, "/maxentResults.csv", sep=""))
      finrow<-length(resultsCV[,1])
      TestAUC<-paste(resultsCV$Test.AUC[finrow])
      Oms<-read.csv(paste(SpOutDirFin, "/",mapsp,"_omission.csv", sep=""))
      Omsn<-which(Oms$Training.omission==max(Oms$Training.omission[Oms$Training.omission<=0.05]))
      myt<-Oms$Corresponding.Cloglog.value[Omsn]
      
      # thresholds
      thresh1<-as.numeric(results[1,c("Maximum.training.sensitivity.plus.specificity.Cloglog.threshold")])
      thresh2<-as.numeric(results[1,c("Equate.entropy.of.thresholded.and.original.distributions.Cloglog.threshold")])
      thresh3<-as.numeric(results[1,c("Balance.training.omission..predicted.area.and.threshold.value.Cloglog.threshold")])
      thresh4<-myt
      
      threshs<-c(thresh1,thresh2,thresh3,thresh4)
      #t<-as.numeric(DatSum$thresh[is])[1] #only use first one as they are all the same for each species within this mapping unit
      t <- 3
      thresh<-threshs[t]
      
      #DatSum$thresh[is]<-thresh #write out threshold used for each species into summary table
      #write.table(x=DatSum, file=paste(MyDir, "/Maxent_LUTables/",projectName,"_Fin3x.csv",sep=""),   
      #            na = "NA", row.names = FALSE, col.names = TRUE, sep=",")
      
      ###########################################
      
      Sys.sleep(0.1)
      print(paste(Sys.time(), "reading in model for modelling unit", mapsp, sep=" "))
      # unzip and read SDM projection = inverted cost raster
      SDMfull<-raster(paste(SpOutFile), crs="+init=epsg:4326")
      SDMfull<-merge(SDMfull,GLOBraster)
      
      
      
      # whole modelling group cost distance
      # ###########################################
      # ###########################################
      # #if this mapspecies includes more than one species, calculate overall cost distance for modeling unit to compare each species to later
       if(length(is)>1){
         
         maprecs<-data.frame()
         for (sp in myspp){
           if(file.exists(paste(OccDir,"/",sp,"_allRecs.csv",sep=""))){
             addrecso<-read.table(paste(OccDir,"/",sp,"_allRecs.csv",sep=""), sep=",", header=TRUE)   #read in species recs
             addrecs<-addrecso[addrecso$uncertainty<= (100000),]
             if(length(addrecs[,1])==0){
               addrecs<-addrecso
             }
             maprecs<-rbind(maprecs,addrecs)
           }}

        if(length(maprecs[,1])>0){

          Sys.sleep(0.1)
          print(paste(Sys.time(), "making extent for modelling unit ",mapsp, sep=" "))

          rangewidthDG<-max((max(na.omit(maprecs$long))-min(na.omit(maprecs$long))),
                            (max(na.omit(maprecs$lat))-min(na.omit(maprecs$lat))))
          rangewidthM<-rangewidthDG*100000 #turn dec deg into meters
          rangewidthM<-ifelse(rangewidthM>3000000,3000000,rangewidthM) #make buffer no bigger than this
          rangewidthM<-ifelse(rangewidthM<1000000,1000000,rangewidthM) #make buffer no smaller than this
          yMOD<-rangewidthM

          #make extent object describing a box around occurrence records (margin= 4 times x and y extent of occurrence records):
          xrange<- (max(na.omit(maprecs$long)))-(min(na.omit(maprecs$long)))
          yrange<- (max(na.omit(maprecs$lat)))-(min(na.omit(maprecs$lat)))
          xcenter<- (min(na.omit(maprecs$long)))+(xrange/2)
          ycenter<- (min(na.omit(maprecs$lat)))+(yrange/2)
          rangemax<-max(xrange,yrange)
          xmin<-ifelse((xcenter-(rangemax/2)- yMOD/100000) >= -180, (xcenter-(rangemax/2)- yMOD/100000), -180)
          xmax<-ifelse((xcenter+(rangemax/2)+ yMOD/100000) <=  180, (xcenter+(rangemax/2)+ yMOD/100000),  180)
          ymin<-ifelse((ycenter-(rangemax/2)- yMOD/100000) >=  -90, (ycenter-(rangemax/2)- yMOD/100000),  -90)
          ymax<-ifelse((ycenter+(rangemax/2)+ yMOD/100000) <=   90, (ycenter+(rangemax/2)+ yMOD/100000),   90)
          myextMOD<-extent(xmin,xmax,ymin, ymax)
        } else{
          myextMOD<-extent(-180,180,-90,90)
        }

  SDMfull <- crop(SDMfull, y=myextMOD)
  AreaRasterB<-crop(GLOBrasterB, myextMOD)
  AreaRaster<-crop(GLOBraster, myextMOD)
  snapraster<-crop(GLOBrasterC, myextMOD)

  ######################################################################################
  # calculate cost distance from currently invaded points = likelihood of occupancy

  Sys.sleep(0.1)
  print(paste(Sys.time(), "calculating cost distance rasters for", mapsp, sep=" "))

  #only continue if there's any points
  coords<-data.frame(long=maprecs$long,lat=maprecs$lat)  # get coordinates of occurrence records
  spcoords <- SpatialPoints(coords[,1:2])              # makes spatial version of occurrence coordinates
  spcoords@proj4string<-CRS("+init=epsg:4326")         # define projection to match my .asc files

  mystack<-stack()

  if(ncell(SDMfull)>130000000){

    myextMOD1<-extent(xmin,xcenter+5,ymin,ycenter+5)
    myextMOD2<-extent(xcenter-5,xmax,ymin,ycenter+5)
    myextMOD3<-extent(xmin,xcenter+5,ycenter-5,ymax)
    myextMOD4<-extent(xcenter-5,xmax,ycenter-5,ymax)
    myextMODs<-c(myextMOD1,myextMOD2,myextMOD3,myextMOD4)

    for (i in c(1:length(myextMODs))){
      thisext<-myextMODs[[i]]
      print(paste("creating partial CD raster for extent",i))
      print(myextMODs[[i]])

      subSDM<-crop(SDMfull,thisext)
      # turn SDM into cost raster (high suitability=low cost to dispersal)
      model_cost.ras  <- -1 * log(subSDM)     # this is the original version of model cost
      model_trans.ras <- 1 / model_cost.ras    # but using the inverse here, to fit the accCost
      # function which is based on a transition matrix. For accCost() the cost of moving between
      # cells = 1 / transition value.

      # create transition matrix based on SDM cost raster:
      gc()
      trans     <- transition(model_trans.ras, transitionFunction=mean, directions = 8)
      trans     <- geoCorrection(trans, type="c", multpl=F)
      gc()

      # calculate accumalted cost & save (= 2.)
      subldr1<-accCost(trans, spcoords)
      subldr2<-merge(snapraster,subldr1)
      # change zero values to a very small non-zero value, to avoid nodata in division
      subldr2[subldr2 < 0.0005 & !(is.na(subldr2)) & subldr2 != -Inf] <- 0.0005
      subldr2[subldr2 == Inf] <-  9999999999

      mystack<-stack(mystack,subldr2)
    }

    lin_dist.rasMU<- calc(mystack,fun=function(x) {min(x,na.rm=TRUE)})

  } else {

    #turn SDM into cost raster (high suitability=low cost to dispersal)
    model_cost.ras  <- -1 * log(SDMfull)     # this is the original version of model cost
    model_trans.ras <- 1 / model_cost.ras    # but using the inverse here, to fit the accCost
    # function which is based on a transition matrix. For accCost() the cost of moving between
    # cells = 1 / transition value.

    # create transition matrix based on SDM cost raster:
    gc()
    trans     <- transition(model_trans.ras, transitionFunction=mean, directions = 8)
    trans     <- geoCorrection(trans, type="c", multpl=F)
    gc()

    # calculate accumalted cost & save (= 2.)
    lin_dist.rasMU = accCost(trans, spcoords)
    # change zero values to a very small non-zero value, to avoid nodata in division
    lin_dist.rasMU[lin_dist.rasMU < 0.0005 & !(is.na(lin_dist.rasMU)) & lin_dist.rasMU != -Inf] <- 0.0005
    lin_dist.rasMU[lin_dist.rasMU == Inf] <-  9999999999

  }

  #save absolute cost distance raster (within x cost distance)
  writeRaster(lin_dist.rasMU, filename=paste(SpOutDirFin,"/ModUnit_",mapsp,"_",cat,"_AbsCostDistance.asc",sep=""), format="ascii", overwrite=TRUE)

 }
# ###########################################

      
   
      ###########################################
      # details for each species in this unit ###
      ###########################################
      
      for (sp in myspp){
        

        ########################################
        #delete any temp files older than 1h
        #removeTmpFiles(h=0.1)
        ########################################
        
        n<-which(DatSum$gen_sp_subtax==sp) #gets row number for current species
        m<-which(myspp==sp)
        print(paste(Sys.time(), "species number", m, "out of", length(myspp), 
                    "(",sp,")","for modelling unit number", mapspnum, 
                    "(",mapsp,")", sep=" "))
        
        #get species records:
        if(file.exists(paste(OccDir,"/",sp,"_allRecs.csv",sep=""))){
          
          sprecso<-read.table(paste(OccDir,"/",sp,"_allRecs.csv",sep=""), sep=",", header=TRUE)   #read in species recs
          sprecs<-sprecso[sprecso$uncertainty<= (100000),]
          if(length(sprecs[,1])==0){
            sprecs<-sprecso
          }
          
          if(length(sprecs[,1])>0){
            
            Sys.sleep(0.1)
            print(paste(Sys.time(), "reading records & making extent for ",sp, sep=" "))
            
            rangewidthDG<-max((max(na.omit(sprecs$long))-min(na.omit(sprecs$long))),
                              (max(na.omit(sprecs$lat))-min(na.omit(sprecs$lat))))
            rangewidthM<-rangewidthDG*100000 #turn dec deg into meters
            actualRWM<-rangewidthM-250000 #make it 250km smaller than extent to avoid harsh straight cut off lines
            actualRWM<-ifelse(actualRWM<250000,250000,actualRWM) #make buffer no smaller than this
            rangewidthM<-ifelse(rangewidthM>3000000,3000000,rangewidthM) #make buffer no biger than this
            rangewidthM<-ifelse(rangewidthM<1000000,1000000,rangewidthM) #make buffer no smaller than this
            y<-rangewidthM
            
            #make extent object describing a box around occurrence records (margin= 4 times x and y extent of occurrence records):
            xrange<- (max(na.omit(sprecs$long)))-(min(na.omit(sprecs$long)))
            yrange<- (max(na.omit(sprecs$lat)))-(min(na.omit(sprecs$lat)))
            xcenter<- (min(na.omit(sprecs$long)))+(xrange/2)
            ycenter<- (min(na.omit(sprecs$lat)))+(yrange/2)
            rangemax<-max(xrange,yrange)
            xmin<-ifelse((xcenter-(rangemax/2)- y/100000) >= -180, (xcenter-(rangemax/2)- y/100000), -180)
            xmax<-ifelse((xcenter+(rangemax/2)+ y/100000) <=  180, (xcenter+(rangemax/2)+ y/100000),  180)
            ymin<-ifelse((ycenter-(rangemax/2)- y/100000) >=  -90, (ycenter-(rangemax/2)- y/100000),  -90)
            ymax<-ifelse((ycenter+(rangemax/2)+ y/100000) <=   90, (ycenter+(rangemax/2)+ y/100000),   90)
            myext<-extent(xmin,xmax,ymin, ymax)
          } else{
            myext<-extent(-180,180,-90,90)
          }
          
          
          # read raw SDM raster and prep reference raster for this area:
          SDMfull<-raster(paste(SpOutFile), crs="+init=epsg:4326")
          SDMfull<-merge(SDMfull,GLOBraster)
          SDM <- crop(SDMfull, y=myext)
          
          AreaRasterB<-crop(GLOBrasterB, myext)
          AreaRaster<-crop(GLOBraster, myext)
          snapraster<-crop(GLOBrasterC, myext)
          
          print(SDM)
          
          ##############################################################################################
          #### cost distance calculations ##############################################################
          ##############################################################################################
          
          ######################################################################################
          # calculate cost distance from currently occupied points
          
          Sys.sleep(0.1)
          print(paste(Sys.time(), "calculating cost distance rasters for", sp, sep=" "))
          
          if(length(sprecs[,1])>0){     #only continue if there's any points
            
            coords<-data.frame(long=sprecs$long,lat=sprecs$lat)  # get coordinates of occurrence records
            spcoords <- SpatialPoints(coords[,1:2])              # makes spatial version of occurrence coordinates
            spcoords@proj4string<-CRS("+init=epsg:4326")         # define projection to match my .asc files
            
            Sys.sleep(0.1)
            print(paste(Sys.time(), "absolute CD", sep=" "))
            
            
            if(ncell(SDM)>130000000){
              
              mystack<-stack()
              
              myext1<-extent(xmin,xcenter+5,ymin,ycenter+5)
              myext2<-extent(xcenter-5,xmax,ymin,ycenter+5)
              myext3<-extent(xmin,xcenter+5,ycenter-5,ymax)
              myext4<-extent(xcenter-5,xmax,ycenter-5,ymax)
              myexts<-c(myext1,myext2,myext3,myext4)
              
              for (i in c(1:length(myexts))){
                thisext<-myexts[[i]]
                print(paste("creating partial CD raster for extent",i))
                print(myexts[[i]])
                
                
                subSDM<-crop(SDM,thisext) 
                # turn SDM into cost raster (high suitability=low cost to dispersal)
                model_cost.ras  <- -1 * log(subSDM)     # this is the original version of model cost
                model_trans.ras <- 1 / model_cost.ras    # but using the inverse here, to fit the accCost
                # function which is based on a transition matrix. For accCost() the cost of moving between
                # cells = 1 / transition value.
                
                # create transition matrix based on SDM cost raster:
                gc()
                trans     <- transition(model_trans.ras, transitionFunction=mean, directions = 8)
                trans     <- geoCorrection(trans, type="c", multpl=F)
                gc()
                
                # calculate accumalted cost & save (= 2.)
                subldr1<-accCost(trans, spcoords)
                subldr2<-merge(snapraster,subldr1)
                # change zero values to a very small non-zero value, to avoid nodata in division
                subldr2[subldr2 < 0.0005 & !(is.na(subldr2)) & subldr2 != -Inf] <- 0.0005
                subldr2[subldr2 == Inf] <-  9999999999
                
                mystack<-stack(mystack,subldr2)
              }
              
              lin_dist.ras<- calc(mystack,fun=function(x) {min(x,na.rm=TRUE)})
              
              
            } else {
              
              # turn SDM into cost raster (high suitability=low cost to dispersal)
              model_cost.ras  <- -1 * log(SDM)     # this is the original version of model cost
              model_trans.ras <- 1 / model_cost.ras    # but using the inverse here, to fit the accCost
              # function which is based on a transition matrix. For accCost() the cost of moving between
              # cells = 1 / transition value.
              
              # create transition matrix based on SDM cost raster:
              gc()
              trans     <- transition(model_trans.ras, transitionFunction=mean, directions = 8)
              trans     <- geoCorrection(trans, type="c", multpl=F)
              gc()
              
              # calculate accumalted cost & save (= 2.)
              lin_dist.ras = accCost(trans, spcoords)
              # change zero values to a very small non-zero value, to avoid nodata in division
              lin_dist.ras[lin_dist.ras < 0.0005 & !(is.na(lin_dist.ras)) & lin_dist.ras != -Inf] <- 0.0005
              lin_dist.ras[lin_dist.ras == Inf] <-  9999999999
              
            }
            
            # ##############################################################################################
            # #### split individual species off from modelling unit ########################################
            # ##############################################################################################
            # 
            # #make any cells that are closer to another species within this modeling unit 9999999999 i.e. really large cost distance
            # if(length(is)>1){
            #   lin_dist.rasMU<-raster(paste(SpOutDirFin,"/ModUnit_",mapsp,"_",cat,"_AbsCostDistance.asc",sep=""), crs="+init=epsg:4326")
            #   thislin_dist.rasMU<-crop(lin_dist.rasMU,extent(lin_dist.ras))
            #   #calculate difference between overall modelling unit cost distance and this species' cost distance -> where values are negative,
            #   #the cost distance for this species is greater than for the modelling unit overall, therefore cost distance to another species'
            #   #record is smaller than for any of this species' records i.e. it belongs to another lineage
            #   rasdif<- thislin_dist.rasMU-lin_dist.ras
            #   values(rasdif)<- ifelse(values(rasdif)< -x1/10,-1,1) #this is a x2/10 km (=50km) CD overlap between distributions
            #   lin_dist.ras<- lin_dist.ras*rasdif
            #   values(lin_dist.ras)<- ifelse(values(lin_dist.ras)< -1, 9999999999,values(lin_dist.ras))
            # }
            # 
            

            #save absolute cost distance raster
            writeRaster(lin_dist.ras, filename=paste(IndSppOutDir,"/",sp,"_",cat,"_AbsCostDistance.asc",sep=""), format="ascii", overwrite=TRUE)
            
            ##############################################################################################
            #### create evolutionarily feasible distance from occurrences (y) for model cutting ##########
            ##############################################################################################
            
            # if(length(is)>1){
            # coords<-data.frame(long=maprecs$long,lat=maprecs$lat)  # get coordinates of occurrence records
            # spcoordsMOD <- SpatialPoints(coords[,1:2],proj4string=CRS("+init=epsg:4326"))
            # Buf<-buffer(spcoordsMOD, width=yMOD, dissolve=TRUE) #make y km buffer around current occurrences
            # }else{
              coords<-data.frame(long=sprecs$long,lat=sprecs$lat)  # get coordinates of occurrence records
              spcoords <- SpatialPoints(coords[,1:2],proj4string=CRS("+init=epsg:4326"))
              Buf<-buffer(spcoords, width=actualRWM, dissolve=TRUE) #make y km buffer around current occurrences
            # }
            
            
            Buf1<-rasterize(x=Buf,y=snapraster,as.numeric(1))
            Buf1<-Buf1*AreaRasterB #keeps anything on land as is (*1) but makes ocean values of buffer NA (*NA)
            Buf1<-merge(Buf1,AreaRaster) #make any non-buffer land 0 again not NA, buffered area stays 1
            
            Sys.sleep(0.1)
            print(paste(Sys.time(), "relative  CD", sep=" "))
            
            ###########################################################################################
            ###### creating relative cost distance rasters ############################################
            ###########################################################################################
            
            # ###################################
            # #X1################################
            # 
            # #convert distance into fraction of max distance and invert to get risk (close=high risk=1; far=low risk=0)
            # #then make anything greater than x1 cost distance km = 0
            # lin_dist.rasX1<-lin_dist.ras/x1
            # lin_dist.rasX1<- 1-lin_dist.rasX1
            # lin_dist.rasX1[lin_dist.rasX1 <= 0 & !(is.na(lin_dist.rasX1))] <- 0 
            # 
            # #get rid of weird huge values (they disappear when cut to current extent so 
            # #must be outside somewhere where it should be 0)
            # values(lin_dist.rasX1)<-ifelse(values(lin_dist.rasX1) != -Inf, values(lin_dist.rasX1), 0)
            # values(lin_dist.rasX1)<-ifelse(values(lin_dist.rasX1) != Inf, values(lin_dist.rasX1), 0)
            # values(lin_dist.rasX1)<- ifelse(values(lin_dist.rasX1)>1,0,values(lin_dist.rasX1))
            # lin_dist.rasX1<-lin_dist.rasX1*AreaRasterB #make oceans NA again
            # 
            # #restrict to y km buffer around occurrences
            # lin_dist.rasX1<-lin_dist.rasX1*Buf1
            # 
            # #save inverted relative cost distance raster (within x1 cost distance)
            # writeRaster(lin_dist.rasX1, filename=paste(IndSppOutDir,"/",sp,"_",cat,"_RelCostDistanceX1.asc",sep=""), format="ascii", overwrite=TRUE)
            # lin_dist.rasCutX1<-lin_dist.rasX1
            # values(lin_dist.rasCutX1)<-ifelse(values(lin_dist.rasCutX1) <= 0 ,0,1)
            # #save binary occurrence raster = anywhere where cost distance <x1 
            # writeRaster(lin_dist.rasCutX1, filename=paste(IndSppOutDir,"/",sp,"_",cat,"_RelCostDistanceBinX1.asc",sep=""), format="ascii", overwrite=TRUE)
            
            ###################################
            #X2################################
            
            #convert distance into fraction of max distance and invert to get risk (close=high risk=1; far=low risk=0)
            #maxy<-max(na.omit(values(lin_dist.ras)[values(lin_dist.ras)!=Inf & values(lin_dist.ras)!= -Inf]))
            #then make anything greater than x cost distance km = 0
            lin_dist.ras<-lin_dist.ras/x2
            lin_dist.ras<- 1-lin_dist.ras
  
            # #new clause, make distance threshold stronger for species that were modeled as group
            # if(length(myspp)>1){
            #   lin_dist.ras[lin_dist.ras < 0.5 & !(is.na(lin_dist.ras))] <- 0 
            # }else{
             lin_dist.ras[lin_dist.ras <= 0 & !(is.na(lin_dist.ras))] <- 0 
            # }
            # 

            #get rid of weird huge values (they disappear when cut to current extent so 
            #must be outside somewhere where it should be 0)
            values(lin_dist.ras)<-ifelse(values(lin_dist.ras) != -Inf, values(lin_dist.ras), 0)
            values(lin_dist.ras)<-ifelse(values(lin_dist.ras) != Inf, values(lin_dist.ras), 0)
            values(lin_dist.ras)<- ifelse(values(lin_dist.ras)>1,0,values(lin_dist.ras))
            lin_dist.ras<-lin_dist.ras*AreaRasterB #make oceans NA again
            
            #restrict to y km buffer around occurrences
            lin_dist.ras<-lin_dist.ras*Buf1
            
            
            
            
            
            
            #save inverted relative cost distance raster (within x2 cost distance)
            writeRaster(lin_dist.ras, filename=paste(IndSppOutDir,"/",sp,"_",cat,"_RelCostDistance.asc",sep=""), format="ascii", overwrite=TRUE)
            lin_dist.rasCut<-lin_dist.ras
            values(lin_dist.rasCut)<-ifelse(values(lin_dist.rasCut) <= 0 ,0,1)
            #save binary occurrence raster = anywhere where cost distance <x2 
            writeRaster(lin_dist.rasCut, filename=paste(IndSppOutDir,"/",sp,"_",cat,"_RelCostDistanceBin.asc",sep=""), format="ascii", overwrite=TRUE)
            
            ###########################################################################################
            ###### create final raster outputs ########################################################
            ###########################################################################################
            
            # Sys.sleep(0.1)
            # print(paste(Sys.time(), "CD weighted SDM", sep=" "))
            # weighSDM<-SDM*lin_dist.ras
            # values(weighSDM)<-ifelse(values(weighSDM)<thresh,0,values(weighSDM)) #threshold 
            # weighSDM<-((weighSDM-thresh)/(1-thresh)) #rescale 0-1
            # writeRaster(weighSDM, filename=paste(IndSppOutDir,"/",sp,"_",cat,"_CropThreshCDWeighted.asc",sep=""), format="ascii", overwrite=TRUE)
            # 
            Sys.sleep(0.1)
            print(paste(Sys.time(), "CD cut SDM", sep=" "))
            
            # #SDMS cut to cost distance x1
            # cutSDM<-SDM*lin_dist.rasCutX1
            # values(cutSDM)<-ifelse(values(cutSDM)<thresh,0,values(cutSDM)) #threshold 
            # cutSDM<-((cutSDM-thresh)/(1-thresh)) #rescale 0-1
            # writeRaster(cutSDM, filename=paste(IndSppOutDir,"/",sp,"_",cat,"_CropThreshCDCutX1.asc",sep=""), format="ascii", overwrite=TRUE)
            # 
            # cutSDMBin<-cutSDM
            # values(cutSDMBin)<-ifelse(values(cutSDMBin)<thresh,0,1) #threshold 
            # writeRaster(cutSDMBin, filename=paste(IndSppOutDir,"/",sp,"_",cat,"_CropThreshCDCutBinX1.asc",sep=""), format="ascii", overwrite=TRUE)
            # 
            
            
            ##########################################################################
            ### if there's several species in this modelling unit, claculate difference in relative distance to this and other species
            ##########################################################################
            
            if(length(is)>1){
              #load modelling unit absolute cost distance:
              lin_dist.rasMU<-raster(paste(SpOutDirFin,"/ModUnit_",mapsp,"_",cat,"_AbsCostDistance.asc",sep=""), crs="+init=epsg:4326")
              lin_dist.rasMU<-crop(lin_dist.rasMU,extent(lin_dist.ras))
              
              #convert into relative cost distance:
              lin_dist.rasMU<-lin_dist.rasMU/x2
              lin_dist.rasMU<- 1-lin_dist.rasMU
              lin_dist.rasMU[lin_dist.rasMU <= 0 & !(is.na(lin_dist.rasMU))] <- 0 

              #get rid of weird huge values 
              values(lin_dist.rasMU)<-ifelse(values(lin_dist.rasMU) != -Inf, values(lin_dist.rasMU), 0)
              values(lin_dist.rasMU)<-ifelse(values(lin_dist.rasMU) != Inf, values(lin_dist.rasMU), 0)
              values(lin_dist.rasMU)<- ifelse(values(lin_dist.rasMU)>1,0,values(lin_dist.rasMU))
              lin_dist.rasMU<-lin_dist.rasMU*AreaRasterB #make oceans NA again
              
              
              
              
              
              #calculate difference in relative cost distance:
              rasdif<- lin_dist.rasMU-lin_dist.ras #areas closer to other species have higher residual values closer to 1
              values(rasdif)<- ifelse(values(rasdif)> 0,(1-values(rasdif))^15,1) #this is a x2/10 km (=50km) CD overlap between distributions from which to fade out
              values(rasdif)<- ifelse(values(rasdif)>1,1,values(rasdif))
              values(rasdif)<- ifelse(values(rasdif)<0,0,values(rasdif))
              SDMfin<-SDM*rasdif
            } else{
              SDMfin<-SDM
            }
            
            #SDMS cut to cost distance x2
            cutSDM<-SDMfin*lin_dist.rasCut
            values(cutSDM)<-ifelse(values(cutSDM)<thresh,0,values(cutSDM)) #threshold 
            cutSDM<-((cutSDM-thresh)/(1-thresh)) #rescale 0-1
            values(cutSDM)<-ifelse(values(cutSDM)<=0,0,values(cutSDM)) #threshold 
            writeRaster(cutSDM, filename=paste(IndSppOutDir,"/",sp,"_",cat,"_CropThreshCDCut.asc",sep=""), format="ascii", overwrite=TRUE)
            
            cutSDMBin<-cutSDM
            values(cutSDMBin)<-ifelse(values(cutSDMBin)<=0,0,1) #threshold 
            writeRaster(cutSDMBin, filename=paste(IndSppOutDir,"/",sp,"_",cat,"_CropThreshCDCutBin.asc",sep=""), format="ascii", overwrite=TRUE)
            
            ########################################################################################
            ##### plotting results #################################################################
            ########################################################################################
            
            #colours:
            #define plot colors
            mycol_terrain = rev(terrain.colors(n = 100, alpha=1))
            mycol_terrainB = rev(terrain.colors(n = 100, alpha=0.5)) #alpha=0:1, higher = less transparent

            #transparency prefixes: 0=00 (transparent), 0.1=1A, 0.2=33, 0.3=4D, 0.4=66, 0.5=80, 0.6=99, 0.7=B3, 0.8=CC, 0.9=E6, 1=FF(solid)
            mycol_blues <-  c(colorRampPalette(c("lemonchiffon","lightsteelblue1","royalblue","navy"))(200))
            mycol_blues <- paste(mycol_blues,"99",sep="")
            mycol_reds <-  c(colorRampPalette(c("lightsteelblue1","lemonchiffon","yellow","orange","orangered","red4"))(200))
            mycol_reds <- paste(mycol_reds,"99",sep="")
            mycol_greys <-  c(colorRampPalette(c("grey95","grey30","grey10","black","black"))(200))
            mycol_greys <- paste(mycol_greys,"99",sep="")
            WHOcol <- rgb(150, 180, 250, max = 255, alpha = 0) #alpha=0 means completely transparent
            mycol_01 <-  c(colorRampPalette(c("white","lemonchiffon"))(200))
            mycol_01 <- paste(mycol_01,"99",sep="")
            
            ###################################
            # prepare data for plots
            
            mybg<-crop(bg, extent(SDM))
            species_shapeA <- subset(WHOpolies, WHOpolies$Scientific == gsub("_"," ",sp))
            
            ######################################################################################
            # make .png plot
            
            #plot:
            
            Sys.sleep(0.1)
            print(paste(Sys.time(),"plotting outputs",sep=" "))
            
            png(filename=paste(PlotDir,"/",sp,"_",cat,"_SDMVet-revised.png",sep=""), units="in", width=10, height=12, res=500) #png smaller file size hat pdf
            par(mfrow = c(3,2), mar = c(1,1,1,1), oma=c(1,1,5,1), mgp=c(1.2,0.5,0))
            #mar=c(bottom, left, top, right); oma= outer margins; mpg = axis label locations c(3, 1, 0)
            
            #1.unthresholded, cropped raw SDM
            
            Sys.sleep(0.1)
            print(paste(Sys.time(),"plot 1",sep=" "))
            
            plot(mybg, col=mycol_greys, ylab="latitude", xlab="longitude", legend=FALSE, cex.axis=0.8, bty="l")
            plot(SDM, col=mycol_blues[c(1,50:200)], legend=TRUE, zlim=c(0,1), add=TRUE,maxpixels=1e8)
            plot(worldmap[5], border="grey20", lwd= 0.5, add=TRUE, col=NA, reset=FALSE)
            
            if(length(species_shapeA)>=1){
              plot(species_shapeA[2], add=TRUE, border="gold", col=WHOcol, lwd= 1)
            }
            
            points(sprecs$lat~sprecs$long, col="black", pch=1, cex=1)
            
            if(length(sprecs[,1])>=1){
              myx<-mean(sprecs$long)
              myy<-min(sprecs$lat)-2
              #graphics::text(x=myx, y=myy, labels = gsub("_"," ",sp), cex = 1, col = mycols2[1])
              graphics::text(x=myx, y=myy, labels = paste("Expert Range Estimate"), cex = 1, col = "gold")
            }

            legend(legend=paste(cat,"Unthresholded Raw Model"),x="top", cex=1, bty="o",bg="white")
            
            #2. cost distance from known occurrences
            
            Sys.sleep(0.1)
            print(paste(Sys.time(),"plot 2",sep=" "))
            
            plot(mybg, col=mycol_greys, ylab="latitude", xlab="longitude", legend=FALSE, cex.axis=0.8, bty="l")
            plot(lin_dist.ras, col=mycol_blues[c(1,50:200)], legend=TRUE, zlim=c(0,1), maxpixels=1e8, add=TRUE)
            plot(worldmap[5], border="grey20", lwd= 0.5, add=TRUE, col=NA, reset=FALSE)
            
            if(length(species_shapeA)>=1){
              plot(species_shapeA[2], add=TRUE, border="gold", col=WHOcol, lwd= 1)
            }
            
            points(sprecs$lat~sprecs$long, col="black", pch=1, cex=1)
            
            legend(legend=paste(cat,"Inverted Relative Cost Distance (RCD) <=",x2/1000,"km"),x="top", cex=1, bty="o",bg="white")
            
            #3. binary SDM cut to areas within x2 maximum cost distance (cut to predefined cost distance)
            
            Sys.sleep(0.1)
            print(paste(Sys.time(),"plot 3",sep=" "))
            
            plot(mybg, col=mycol_greys, ylab="latitude", xlab="longitude", cex.axis=0.8, legend=FALSE, bty="l")
            plot(cutSDMBin, col=mycol_blues[c(1,150)], legend=FALSE, add=TRUE,maxpixels=1e8)
            plot(worldmap[5], border="grey20", lwd= 0.5, add=TRUE,col=NA,reset=FALSE)
            
            if(length(species_shapeA)>=1){
              plot(species_shapeA[2], add=TRUE, border="gold", col=WHOcol, lwd= 1)
            }
            
            legend(legend=paste(cat,"Binary Model Cut to RCD<=", x2/1000, "& Thresholded"),x="top", cex=1, bty="o",bg="white")
            
            #4. SDM cut to areas within x2 maximum cost distance (cut to predefined cost distance)
            
            Sys.sleep(0.1)
            print(paste(Sys.time(),"plot 4",sep=" "))
            
            plot(mybg, col=mycol_greys, ylab="latitude", xlab="longitude",cex.axis=0.8, legend=FALSE, bty="l")
            plot(cutSDM, col=mycol_blues[c(1,50:200)], legend=TRUE, zlim=c(0,1), add=TRUE, maxpixels=1e8)
            plot(worldmap[5], border="grey20", lwd= 0.5, add=TRUE,col=NA,reset=FALSE)
            
            if(length(species_shapeA)>=1){
              plot(species_shapeA[2], add=TRUE, border="gold", col=WHOcol, lwd= 1)
            }
            
            legend(legend=paste(cat,"Continuous Model Cut to RCD<=", x2/1000, "& Thresholded"),x="top", cex=1, bty="o",bg="white")
            
            #5. current SDM + Buffer
            
            Sys.sleep(0.1)
            print(paste(Sys.time(),"plot 5",sep=" "))
            
            lin_dist.ras<-lin_dist.ras*cutSDMBin
            #values(lin_dist.ras)<-ifelse(values(lin_dist.ras)<=0,NA,values(lin_dist.ras))
            values(lin_dist.ras)<-ifelse(values(lin_dist.ras)>0 & values(lin_dist.ras) <0.01,0.011,values(lin_dist.ras))
            myvals<-values(lin_dist.ras)
            myvals<-myvals[myvals>0]
            CDQ75<-round(quantile(myvals,0.75, na.rm=TRUE),2)
            CDQ75<-ifelse(CDQ75==1,0.99,CDQ75)
            CDQ75<-ifelse(CDQ75< 0.03,0.03,CDQ75)
            print(paste("CDQ75 =",CDQ75))
            CDQ50<-round(quantile(myvals,0.5, na.rm=TRUE),2)
            CDQ50<-ifelse(CDQ50==1,CDQ50-0.02,CDQ50)
            CDQ50<-ifelse(CDQ50==CDQ75,CDQ50-0.01,CDQ50)
            CDQ50<-ifelse(CDQ50< 0.02,0.02,CDQ50)
            print(paste("CDQ50 =",CDQ50))
            CDQ25<-round(quantile(myvals,0.25, na.rm=TRUE),2)
            CDQ25<-ifelse(CDQ25==1,CDQ25-0.03,CDQ25)
            CDQ25<-ifelse(CDQ25==CDQ75,CDQ25-0.02,CDQ25)
            CDQ25<-ifelse(CDQ25==CDQ50,CDQ25-0.01,CDQ25)
            CDQ25<-ifelse(CDQ25<0.01,0.01,CDQ25)
            print(paste("CDQ25 =",CDQ25))
            if(is.na(CDQ75)){
              CDQ75<-0.75
              CDQ50<-0.5
              CDQ25<-0.25
              print(paste("default CDQ breaks = 0, 0.25, 0.5, 0.75, 1"))
            }
            
            plot(mybg, col=mycol_greys, ylab="latitude", xlab="longitude",cex.axis=0.8, legend=FALSE, bty="l")
            #plot(Buf1, col=mycol_01[c(1,200)], legend=FALSE, maxpixels=1e8, colNA="lightgrey",add=TRUE)
            plot(lin_dist.ras, col=mycol_blues[c(1,50,100,150,200)], legend=TRUE, zlim=c(0,1), add=TRUE, maxpixels=1e8, breaks=c(0,0.0001,CDQ25,CDQ50,CDQ75,1))
            plot(worldmap[5], border="grey20", lwd= 0.5, add=TRUE, col=NA, reset=FALSE)
            
            if(length(species_shapeA)>=1){
              plot(species_shapeA[2], add=TRUE, border="gold", col=WHOcol, lwd= 1)
            }
            
            legend(legend=paste("Quartiles of Cost Distance (0= Far, 1= Close)"), x="top", cex=1, bty="o", bg="white")
            
            #6. final SDM summary
            
            Sys.sleep(0.1)
            print(paste(Sys.time(),"plot 6",sep=" "))
            
            #values(cutSDM)<-ifelse(values(cutSDM)<=0,NA,values(cutSDM))
            values(cutSDM)<-ifelse(values(cutSDM)>0 & values(cutSDM) <0.01,0.011,values(cutSDM))
            myvals<-values(cutSDM)
            myvals<-myvals[myvals>0]
            SDMQ75<-round(quantile(myvals,0.75, na.rm=TRUE),2)
            SDMQ75<-ifelse(SDMQ75==1,0.99,SDMQ75)
            SDMQ75<-ifelse(SDMQ75< 0.03,0.03,SDMQ75)
            print(paste("SDMQ75 =",SDMQ75))
            SDMQ50<-round(quantile(myvals,0.5, na.rm=TRUE),2)
            SDMQ50<-ifelse(SDMQ50==1,SDMQ50-0.02,SDMQ50)
            SDMQ50<-ifelse(SDMQ50==SDMQ75,SDMQ50-0.01,SDMQ50)
            SDMQ50<-ifelse(SDMQ50< 0.02,0.02,SDMQ50)
            print(paste("SDMQ50 =",SDMQ50))
            SDMQ25<-round(quantile(myvals,0.25, na.rm=TRUE),2)
            SDMQ25<-ifelse(SDMQ25==1,SDMQ25-0.03,SDMQ25)
            SDMQ25<-ifelse(SDMQ25==SDMQ75,SDMQ25-0.02,SDMQ25)
            SDMQ25<-ifelse(SDMQ25==SDMQ50,SDMQ25-0.01,SDMQ25)
            SDMQ25<-ifelse(SDMQ25<0.01,0.01,SDMQ25)
            print(paste("SDMQ25 =",SDMQ25))
            if(is.na(SDMQ75)){
              SDMQ75<-0.75
              SDMQ50<-0.5
              SDMQ25<-0.25
              print(paste("default SDMQ = 0.25, 0.5, 0.75"))
              
            }
            
            plot(mybg, col=mycol_greys, ylab="latitude", xlab="longitude",cex.axis=0.8, legend=FALSE, bty="l")
            #plot(Buf1, col=mycol_01[c(1,200)], legend=FALSE, maxpixels=1e8, colNA="lightgrey",add=TRUE)
            plot(cutSDM, col=mycol_blues[c(1,50,100,150,200)], legend=TRUE, add=TRUE, zlim=c(0,1), maxpixels=1e8, breaks=c(0,0.0001,SDMQ25,SDMQ50,SDMQ75,1))
            plot(worldmap[5], border="grey20", lwd= 0.5, add=TRUE,col=NA,reset=FALSE)
            
            #plot WHO ranges for species in this modelling unit
            
            if(length(species_shapeA)>=1){
              plot(species_shapeA[2], add=TRUE, border="gold", col=WHOcol, lwd= 1)
            }
            
            
            legend(legend=paste("Quartiles of Suitable Cells (0= Unsuitable, 1= Very suitable)"),x="top", cex=1, bty="o",bg="white")
            ListSpecies<-DatSum$Listing_Name[n]
            #plot title in top margin of page:
            mtext(paste(gsub("_"," ",sp),sep=""), font= 4, side = 3, line = 1, outer = TRUE, cex= 2)
            #mtext(paste("(WHO listed as '",ListSpecies,"')",sep=""), font= 3, side = 3, line = -0.4, outer = TRUE, cex= 1)
            
            if(length(myspp)>1){
              if(length(myspp)<=3){
                mtext(paste("(modelled as ", paste(gsub("_"," ",DatSum$ModUnit[n])),", includes: ", paste(gsub("_"," ",c(myspp[1:3])), collapse = ', '),", etc.)",sep=""), font= 3, side = 3, line = -0.4, outer = TRUE, cex= 1)
              }else{
                mtext(paste("(modelled as ", paste(gsub("_"," ",DatSum$ModUnit[n])),", includes: ", paste(gsub("_"," ",c(myspp[1:3])), collapse = ', '),", etc.)",sep=""), font= 3, side = 3, line = -0.4, outer = TRUE, cex= 1)
              }}
            
            
            dev.off()
            
          }}}}
    
    print(paste(Sys.time()," zipping file ",SpOutFile, sep=""))
    
    
    if(file.exists(SpOutFile) & file.exists(paste(SpOutFile,".gz",sep=""))){
      file.remove(SpOutFile)
    } 
    
    if(file.exists(SpOutFile) & !file.exists(paste(SpOutFile,".gz",sep=""))){
    gzip(filename=paste(SpOutFile),overwrite=TRUE, remove=TRUE)
    }

  }}}

# ######################################################
# #delete the temp dir for this scriptL
# unlink(paste(MyDir,"/tmpCutCD",sep=""), recursive = TRUE) 
# ######################################################

############################################
############################################
# next: 9-VAP_Hotspots #####################
############################################
